DEPENDENCE OF MAPK MEDIATED SIGNALING ON ERK ISOFORMS 
AND DIFFERENCES IN NUCLEAR SHUTTLING 

HEATHER A. HARRINGTON*\ MICHAL KOMOROWSKI\ MARIANO BEGUERISSE DIAZ 2, 
GIAN MICHELE RATTO ^ AND MICHAEL P.H. STUMPF^ 



Abstract. The mitogen activated protein kinase (MAPK) family of proteins is involved in regu- 
lating cellular fate activities such as proliferation, differentiation and apoptosis. Their fundamental 
importance has attracted considerable attention on different aspects of the MAPK signaling dy- 
namics; this is particularly true for the Erk/Mek system, which has become the canonical example 
for MAPK signaling systems. Erk exists in many different isoforms, of which the most widely 
studied are Erkl and Erk2. Until recently, these two kinases were considered equivalent as they 
differ only subtly at the sequence level; however, these isoforms exhibit radically different traffick- 
ing between cytoplasm and nucleus. Here we use spatially resolved data on Erkl/2 to develop 
and analyze spatio-temporal models of these cascades; and we discuss how sensitivity analysis can 
be used to discriminate between mechanisms. We are especially interested in understanding why 
two such similar proteins should co-exist in the same organism, as their functional roles appear to 
be different. Our models elucidate some of the factors governing the interplay between processes 
and the Erkl/2 localization in different cellular compartments, including competition between iso- 
forms. This methodology is applicable to a wide range of systems, such as activation cascades, 
where translocation of species occurs via signal pathways. Furthermore, our work may motivate 
additional emphasis for considering potentially different roles for isoforms that differ subtly at the 
sequence level. 
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Introduction 

Cellular decision making processes require balanced and nuanced responses to environmental, 
physiological and developmental signals. Especially in eukaryotes this involves the interplay and 
concerted action of a number of molecular players, that receive external signals, broadcast them 
further into the cytoplasm, and, if a transcriptional response is called for, shuttle into the nucleus 
and activate the transcriptional machinery. The mitogen activated protein kinases (MAPK) are 
important relays in many of these signal transduction processes. They are, for example, involved 
in regulating cellular fates such as proliferation, differentiation and apoptosis [Ij. The most widely 
studied MAPK, Erkl/2, is activated through phosphorylation by Mek, its MAPK kinase (MAPKK), 
and Mek in turn is activated by its cognate kinase; ultimately, the activator of the cascade of MAPK 
kinases is an initiator signal, such as Ras (MAPKK kinase). Activated kinases, such Erk or Mek, 
are deactivated via dephosphorylation by their respective phosphatases. 

In the active state Erk can shuttle into the nucleus where it can induce the relevant cellular 
responses. But continuous nucleo-cytoplasmic shuttling of Mek/Erk has been demonstrated even 
in the absence of activation [21 |3l HI [5]. Mek is believed to play a central role in the shuttling 
mechanism of Erk using its nuclear export sequence; however, due to its nuclear exclusion sequence, 
Mek is found mostly in the cytoplasm whereas Erk can accumulate in the nucleus 1^ . 
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The two distinct isoforms of Erk, Erkl and Erk2, appear to have distinctly different biologi- 
cal characteristics. Transgenic gene knockout mice lacking Erk2 (Erk2~/~) result in embryonic 
lethality, whereas Erkl~/~ are viable, fecund and suffer from only subtle deficiencies, e.g., reduced 
T cell development [H [?]. Erk2 can compensate in Erkl deficient mice through as yet unclear 
mechanisms. Although Erkl/2 share approximately 85% of the same amino acid sequence, both 
activated by Mek, and considered to behave similarly, shuttling of Erkl from cytoplasm to nucleus 
is approximately three times slower than that of Erk2, regardless of whether the cell is starved 
or stimulated [51 [5]. The rate of shuttling results from differences in the N-terminal domain of 
Erkl and Erk2 [9]. In addition to different trafficking speeds, Erkl is approximately four times less 
abundant than Erk2 in NIH 3T3 cells [lOj. It is also believed that competition exists between Erk 
isoforms in binding to Mek |11) . 

Activation cascades have been studied extensively, but given that experimental results often 
rely on incomplete or indirect observations of the underlying dynamics, and in light of sometimes 
partially contradictory experimental results, mathematical models have become increasingly impor- 
tant for understanding these cascades. Published models of activation cascades, such as MAPK, 
typically use a three tier structure, where the top tier is activated by a single phosphorylation, 
and the remaining two tiers are activated via double phosphorylation cycles. This structure can 
incorporate positive feedback loops and the most widely studied model by Huang and Ferrell pre- 
dicted ultrasensitivity in the MAPK signal [12^ [T3t [Hj . Their work initiated the development of 
several alternative MAPK cascade models, which exhibit cell switching [15| [TBI IT7] and provide 
conditions for bistability (e.g., via double phosphorylation through a distributive mechanism [18] 
or MAPKK sequestration [l9]). Subsequent models have also included inhibition by the doubly 
phosphorylated MAPK onto the initial signal Ras, either directly or implicitly between cycles, via 
a negative-feedback loop. These models predicted sustained oscillations [201 [21], which were then 
verified by further modeling and experimental work [221 . Additionally, other studies have found 
that parameter regions affect the dynamical behavior of the system which may yield multi-stability 
or oscillations [23[ [^ [23] . Although many models resemble the form of the Huang-Ferrell model, 
many authors have explored the cascade in a linear manner |2B1 [?71 [2H] in order to study the effects 
of cascade length, signal transduction and adaptation, in the presence of MAPK specific feedback 
structures. 

Shuttling may play a role in the activation/deactivation activity of Mek/Erk in both the nucleus 
and cytoplasm, which has motivated the recent development of spatio-temporal models [291 [21] 
[22l l30l [3T] . A compartment model of MAPK using single phosphorylation and dephosphorylation 
loops [29] was extended to include double phosphorylation of kinases and phosphatases which 
dephosphorylate the kinases [22] . Models have also been developed to study the differences between 
Erk isoforms [321 [5] [23]. A tiered MAPK model with Erkl, Erk2, Mek and upstream species of the 
MAPK cascade was developed identifying that Erk isoforms cross-regulate each other and may be 
linked to cell fate decisions [32]. Another model focused on Erkl dimerization using data from Erkl 
wild type and the Erkl-A4 mutant (unable to dimerize) and improved the nuclear transport model 
of [29[ by exploring different methods of phosphorylation [24] . Finally, a simple model including 
only Erkl and Erk2 species was constructed to study the nucleo-cytoplasmic shuttling speeds |5]. 

The existence of two genetically nearly identical isoforms of Erk is hard to explain, especially 
in light of the above results. Here we seek to understand potential dynamical reasons for the 
existence of the Erkl and Erk2 isoforms. To do so we have to consider their interactions with Mek 
and nuclear/cytoplasmic shuttling, which appear to be important in Mek/Erk mediated cellular 
decision making processes. We integrate these different processes into a single model and compare it 
to recent spatio-temporally resolved data and consider and contrast the behavior in two cell types, 
HeLa and NIH 3T3, and under different conditions. We find that shuttling, i.e., hitherto often 
ignored spatial aspects, exert great influence over the cellular response phenotypes and are able to 



2 



provide a rationale for the existence of different Erk isoforms. Furthermore, our model supports 
the recently posed hypothesis [TT] of the importance of competition between the Erk isoforms. 

This study also lays the ground work for understanding a complex system in the biological 
sciences by developing an appropriate abstractions; in this case, through the development and 
analysis of a minimal model. In addition, these approaches are suitable in contexts where transport 
and/or different roles of molecules are unclear. 

Materials and Methods 

A spatio-temporal model of Erk signalling. A reaction network of the spatio-temporal model 
is built around three species: Erkl, Erk2, and Mek are allowed to exist in different activation states, 
as complexes and in different cellular compartments (see Fig. lA and find equations in Supplemen- 
tary Material). We assume Erk isoforms (Erkl/Erk2) exist in one of three states: inactive (Erk), 
bound to Mek (M-Erk), or active (Erkp). Erk transitions from inactive to active through Mek 
whereby free phosphorylated Mek (Mek) reversibly interacts with inactive Erk. Upon Mek-Erk 
binding and formation of an intermediary complex (M-Erk), Mek phosphorylates Erk resulting 
in activated Erk (Erkp) and then dissociates from Erkp. In order for Erkp to revert back to its 
inactive form, it must undergo dephosphorylation by a phosphatase. To denote the subcellular 
compartment of Erk, species are denoted by a subscript n to indicate nuclear localization, i.e., 
Erkln is Erkl in the nuclear compartment. The biochemical reactions of the full system are given 
in Table 1 with a description in Table 3. Shuttling rates of Erk in HeLa cells have previously been 
established; however, to determine the shuttling of Erk species in NIH 3T3 cells, we performed 
FRAP experiments, as shown in Fig. 2, and these data enabled us to estimate transport rates. 

Cell culture and transfection. NIH 3T3 cells were plated on glass disks cultured in modified 
Dulbecco's medium supplemented with 10% fetal bovine serum, and were transfected using Lipo- 
fectamine (GenePORTER 2, Genlantis, San Diego, CA). To inactivate the ERK pathway, cells were 
starved by keeping them in 1% FBS for 24 h before imaging. Further details and the design of the 
ERK-GFP fusion proteins are found in [5l[8]. 

FRAP measure of nucleus-cytoplasm shuttling. Photo bleaching was preceded by the acqui- 
sition of a pre-bleach image that was used to estimate the loss of fluorescence due to bleaching and 
for data normalization. The nucleus of the cell was photobleached by repeated scans of the nucleus 
at high laser power. Bleaching was applied for approximately 8 s, which was sufficient to quench 
most of the nuclear fluorescence. Bleaching was followed by time-lapse acquisition to measure the 
recovery. Fluorescence was normalized by 

F(t) = ^Nucif) 

where F^ot indicates the fluorescence (corrected for background) measured before bleaching in the 
entire cell and -Fnuc indicates the same measure in the nucleus. This normalization corrects for 
bleaching caused by imaging. All imaging experiments were performed on a Olympus Fluoview 300 
and Leica SL confocal microscopes. 

Model parameterization. We focus our analysis on HeLA cells and mouse NIH 3T3 cells. Phos- 
phorylation, shuttling rates and initial conditions have previously been specified for HeLa cells 
|29^ [33] whereas shuttling rates here have been measured in mouse NIH 3T3 cells O [8] . We use the 
kinetic phosphorylation and dephosphorylation parameters for an initial estimate for Erk2 from the 
compartment model of Erk2 in HeLa cells using a single phosphorylation loop [29]. Erkl phospho- 
rylation was assumed to occur slightly faster than Erk2 in primary erytheroid progenitor (CFU-E) 
cells |32) and we adjusted the baseline parameter values for the complete model accordingly (Ta- 
ble 1). 
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Figure 1. Model construction, (a) Schematic of full Erk/Mek shuttling model. Inactive Erk 
can bind to Mek to form the M-Erk complex. Mek phosphorylates Erk into Erkp, then Mek 
dissociates and Erkp can be dephosphorylated to Erk. This behavior is symmetric for both Erkl 
and Erk2 isoforms. These reactions (solid arrows) can occur in either the cytoplasm (yellow region) 
or nucleus (orange region) and all model species (gray boxes) can shuttle between compartments 
(dashed arrows), (b) Time course of nuclear Erklpn and Erk2pn that reaches steady-state for 
HeLa and NIH 3T3 cells, (i) Time course at baseline conditions, (ii) Time course at Limited Mek 
(Mek = 0.2). (iii) Time course at equal Erkl/2 (Erkl = Erk2 = 0.5). 



Both HeLa and NIH 3T3 cells are 15-20 microns in diameter and have approximately the same 
volume. The import and export shuttling rates of species between the cytoplasm and nucleus include 
the ratio of cytoplasmic to nuclear volume as used in other compartment models [291 [221 El] • We 
use the Erk2 shuttling rates in HeLa cells from [29] and then calibrated the shuttling rates to satisfy 
that Erkl shuttles three times slower than Erk2 as seen in Fig. 2 and [5]. In starved cells there is 
no phosphorylation by Mek and primarily inactive Erk is shuttled; stimulated cells, on the other 
hand, (which have Mek/Erk activation) shuttle predominantly phosphorylated (Erkp) species. The 
shuttling time constants (r) were calculated by first photobleaching the nucleus and then measuring 
the recovery of fluorescence of Erk in NIH 3T3 cells (see [5l [8] for more details and Supplementary 
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Figure 2. Measure of the nucleus-cytoplasm trafficking of Erk-GFP in living cells, a) Imaging of 
NIH 3T3 fibroblasts transfected with Erk2-GFP. Trafficking rates have been measured by analyzing 
the speed of recovery of fluorescence after bleaching of the nucleus (red dots in the first image) : the 
faster is the recovery the higher is the rate of trafficking. Nuclear fluorescence was bleached at time 
and the time lapse sequence shows its recovery. Cells were starved in 1% fetal bovine serum for 24 
hrs before performing the experiment. After compete recovery the cell were stimulated with 10% 
serum and, after allowing 15 min for maximal Erk activation, a second bleach-recovery sequence 
was acquired. The calibration bar is 10 microns, b) Time course of the fluorescence recovery 
acquired for the illustrated cell before (red) and after activation of the Erk pathway (green) . The 
continuous lines are fitting exponentials of time constants 152 and 49 s respectively, c) Erkl and 
Erk2 have different trafficking rates: The four exponentials indicate the normalized time course of 
the fluorescence recovery for Erkl (dotted lines) and Erk2 (continuous lines) in the starved condition 
(red) and after stimulation (green). Rates have been obtained from averaging rates measured in 
30-60 cells in each condition (5). 



Material for transport estimation) . Transport reactions and shuttling rates can be found in Table 2 
and a description in Table 3. 

The non-zero baseline initial concentrations are the cytoplasmic species Erkl, Erk2, and Mek. 
Estimates of initial levels of Erk in previous MAPK models of HeLa cells ranged between 0.96 — 2.1 
[ 33\ [29] and most previous models have initial Mek and Erk concentrations varying between 
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Table 1. Reactions of Erkl/2 and Mek in full model. Each reaction highlights whether it is a 
forward or reversible reaction by the arrows as illustrated by solid arrows in Fig. lA. All reaction 
rates are from |29) . 



Number 


Reaction 


Forward rate (/iM ^ s ^) 


Reverse rate (s i) 


1 


Erkl + Mek ^ M-Erkl 


8.8 X 10-1 


8.8 X 10-2 


2 


M-Erkl Mek + Erklp 


3 X IQ-^ 




3 


Erkl„ + Mek„ ^ M-Erkl„ 


8.8 X 10-1 


8.8 X 10-2 


4 


M-Erkln -> Mekn + Erklpn 


3 X 10-1 




5 


Erk2 + Mek ^ M-Erk2 


8.8 X 10-1 


8.8 X 10-2 


6 


M-Erk2 -> Mek + Erk2p 


2 X 10-1 




7 


Erk2„ + Mekn ^ M-Erk2„ 


8.8 X 10-1 


8.8 X 10-2 


8 


M-Erk2„ Mekn + Erk2pn 


2 X 10-1 




9 


Erklp Erkl 


1.4 X 10-2 




10 


Erklpn Erkln 


1.4 X 10-2 




11 


Erk2p Erk2 


1.4 X 10-2 




12 


Erk2p„ Erk2n 


1.4 X 10-2 





Table 2. Shuttling rates for HeLa and NTH 3T3 cells. Shuttling rates for HeLa and NIH 3T3 cells 
as shown dotted arrows in Fig. 1 A. Parameters are estimated from (2^ and * denotes rate from [5] . 







Baseline HeLa cells 


Baseline NIH 3T3 cells 


Number 


Transport 


Export Rate (s i) 


Import Rate (s i) 


Export Rate (s i) 


Import Rate (s i) 


1 


Erkl ^ Erkl„ 


1.3 X 10-2 


4 X 10"^ * 


2.1 X IQ-^ * 


1.4 X 10-3 * 


2 


Mek ^ Mek„ 


5.4 X 10-1 


4 X 10-2 


5.4 X 10-1 


4 X 10-2 


3 


M-Erkl ^ M-Erkl„ 


2.6 X 10-1 


1.17 X 10-2 * 


2.6 X 10-1 


1.17 X 10-2 * 


4 


Erklp ^ Erklpn 


1.3 X 10-2 


4 X 10-3 * 


5.2 X 10-3 * 


4.7 X 10-3 * 


5 


Erk2 ^ Erk2„ 


1.8 X 10-2 


1.2 X 10-2 


7.8 X 10-3 * 


5.2 X 10-3 * 


6 


M-Erk2 ^ M-Erk2„ 


2.6 X 10-1 


3.5 X 10-2 


2.6 X 10-1 


3.5 X 10-2 


7 


Erk2p ^ Erk2p„ 


1.3 X 10-2 


1.1 X 10-2 


1.65 X 10-2 * 


1.5 X 10-2 * 



0.1 — 3 ijM. ^29j. In HeLa cells, the ratio of Mek to Erk at the initial conditions is approximately 
1.46 |33ll29j . We do not include inactive Mek in this model and we estimate that the ratio of active 
phosphorylated Mek to total Erk is one. Since Erk in previous models had been assumed to be 
approximately 1 //M [33', we assume the relation total (phosphorylated) Mek ~ total Erkl + 
total Erk2, and the baseline initial concentrations of Erk also satisfy the condition that Erkl is 
four times less abundant than Erk2 in NIH 3T3 cells [TU]. We also take into account that HeLa 
cells may have a different ratio of Erkl /Erk2 or different concentrations of active Mek either due to 
initial activation of the MAPK cascade or Mek sequestration effects, etc., than the baseline initial 
conditions as given in Table 4. 

Steady-state analysis. We investigate the behavior of the system for the relevant baseline pa- 
rameters and find that only one permissible and biologically relevant steady-state exists. Further 
exploration of the parameter space of the models using latin-hypercube sampling [3lj confirms that 
there is only one biologically permissible steady-state across the biophysically plausible parameter 
space. This finding of the minimal model is in line with previous studies that looked into the 
classification of single and double phosphorylation cycle networks in terms of whether they yield 
multistability [55] . 

Sensitivity analysis. To study the effects of the parameter values on the dynamics of Erklpn 
and Erk2pn) we performed a detailed and comprehensive sensitivity analysis. We used a standard 
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Table 3. Description of model parameters. 



Parameters 


Description 


fci 


forward interaction rate between cytoplasmic Mek and Erkl 




reverse interaction rate between cytoplasmic Mek and Erkl 


k2 


phosphorylation rate of cytoplasmic Erkl 




forward interaction rate between nuclear Mek and Erkl 


k-3 


reverse interaction rate between nuclear Mek and Erkl 




phosphorylation rate of nuclear Erkl 


h 


forward interaction rate between cytoplasmic Mek and Erk2 




reverse interaction rate between cytoplasmic Mek and Erk2 




phosphorylation rate of cytoplasmic Erk2 


kj 


forward interaction rate between cytoplasmic Mek and Erk2 




reverse interaction rate between cytoplasmic Mek and Erk2 


h 


phosphorylation rate of nuclear Erk2 


kg 


dephosphorylation rate of cytoplamsic Erklp 


kio 


dephosphorylation rate of nuclear Erklp 


fell 


dephosphorylation rate of cytoplasmic Erk2p 


ki2 


dephosphorylation rate of nuclear Erk2p 


Sl 


import of Erkl into the nucleus 


S-l 


export of Erkl out of the nucleus 


32 


import of Mek into the nucleus 


S2 


export of Mek out the nucleus 


S3 


import of M-Erkl into the nucleus 


S^3 


export of M-Erkl into the nucleus 


S4 


import of Erklp into the nucleus 


S4 


export of Erklp out of the nucleus 


S5 


import of Erk2 into the nucleus 


S-5 


export of Erk2 out of the nucleus 


S6 


import of M-Erk2 into the nucleus 


S-6 


export of M-Erk2 out of the nucleus 


S7 


import of Erk2p into the nucleus 


S-T 


export of Erk2p out of the nucleus 



Table 4. Initial conditions of the model variables. Initial conditions of the model variables for 
different cell types. Some species initial conditions may differ due to cell variability or cell type. 
Baseline values estimated from [291 [3]. 



Initial concentration (/iM) 
Species Baseline Equal Erk 1/2 Limited Mek 

(L2 05 02 

Erk2 0.8 0.5 0.8 

Mek 1 1 0.2 



method based on the Fisher information [SB] , which assesses the change in system output as param- 
eters are varied. Suppose that we are interested in the values Erklpn(i) and Erk2pn(t), denoting 
concentrations of nuclear Erklp and Erk2p respectively, at times ti,...,tN. Then the sensitivity 
coefficient for a parameter ki is 

^ /(9Erklpn(tj)\^ fdErk2p^{tj 



(1) s,. = Y.( — + 



dkj J \ dti 



As the trajectories x.{t) are given by solutions of ordinary differential equations their derivatives 
^Q]^^^ can be easily calculated [37] . We have used this type of sensitivity coefficients as they have a 
clear information geometric interpretation as the infinitesimal distance in the state space, behave 
locally as the Kullback-Leibler divergence [38], and allow us to determine directions of trajectories 
change resulting from parameters changes. Treating the initial conditions as parameters provides 
us with their sensitivities in an equally simple way [36]. 



Assuming a constant stimulus, as is the case for our experimental set-up, we simulate the model 
equations (see (Tables 1-4) for details) until steady-state. The time course of active nuclear Erk is 
given in Fig. la, and total nuclear Erkl and total nuclear Erk2 abundances are shown in Fig. Sla 
in the Supplementary Material. We take into account cell variability by studying the behavior 
at different initial concentrations of Erk and Mek within the known constraints; activation of the 
system is determined by the initial condition of active Mek. Mek then activates Erk, which can 
then traffic in and out of the nucleus; we are therefore interested in both the transient activation 
time (time until a species reaches steady-state) and the realised steady-state values of Erklpn and 
Erk2pn (phopho- nuclear Erk isoforms). We consider the realistic scenario of a smaller stimulus 
activating the system and take this into account by decreasing the initial level of activated Mek, 
which decreases the steady-state values of Erkp (Fig. lb)). 

The overarching question is really why do nearly identical Erk 1/2 isoforms exist? In particular we 
seek to address this question by considering two related questions in more detail: (1) How does the 
behavior of these isoforms differ in different cell-types and under different physiological conditions? 
(2) How does nuclear shuttling differ between the two isoforms and how does this affect the overall 
nuclear activity levels of Erk? More generally, we develop a flexible mechanistic description of 
Erk mediated signaling that includes the nuclear translocation dynamics. This framework can 
be adapted to other molecules, e.g., NFAT, Stat3/5 ^39j etc where we have distinct biochemical 
isoforms with different nuclear shuttling dynamics. 

Minimal model. To address the above questions, the behavior of Erk isoforms under vary- 
ing physical conditions is more closely studied in the absence of transport (Fig. 3a). Sensi- 
tivity analyses identify Erklp and Erk2p concentrations as most sensitive to the phosphoryla- 
tion/dephosphorylation (/cg/fen) rates of Erkl and Erk2 (Fig. 3b). Analysis of the initial conditions 
also show that the available phosphorylated Mek increases in importance at lower signal strength 
(limited Mek) (Fig. 3b, c), which in an experimental context, is what we expect. The percentage 
of Mek that is usually phosphorylated depends on signal strength and this enzyme would greatly 
affect the steady state of Erkp. The balance between phosphorylation/dephosphorylation rates and 
initial conditions provides a setting for a competition scheme between Erk isoforms for its activator 
Mek; furthermore, available evidence shows that there is a close match between Mek and Erkl/2 
concentrations. 

We investigate the effects of the phosphorylation/dephosphorylation rates on the relative abun- 
dances of Erk2p to Erklp at baseline Mek and limited Mek conditions. By defining a log-scale 
indicator 



the sign of (p signifies the dominance of a specific phopho-Erk isoforms as physical conditions are 
changed. For example, given that the value of ip identifies the change of phospho-Erk isoforms 
by orders of magnitude, a value of one indicates ten times more Erk2|3 than Erklp whereas a 
negative value of —1 establishes there is ten times more Erklp than Erk2p at steady state. The 
phosphorylation rates and dephosphorylation rates are each varied separately from 10~^ to 10^ in 



Results 
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Figure 3. Minimal Model, (a) Schematic of biochemical Erk/Mek model, (b) Sensitivity analysis 
for the biochemical model. Mean sensitivity index given by violet line (see Materials and methods), 
(i) Sensitivity to parameters, (ii) Sensitivities to initial concentration. 



line with biophysical considerations and is evaluated (Fig. 4a,b). As the rate of phosphorylation 
of Erkl increases, values of (p become negative indicating dominance by Erklp (see bottom right 
corner of Figs. 4a (ii)), whereas an increase in the Erk2 phosphorylation rate demonstrates Erk2 
prevalence. The effect of the phosphorylation rate on is asymmetric, meaning that there is a larger 
parameter space of values giving rise to Erk2p dominance (see non-overlapping curves in Fig. 4a 
(i)). Notably, the competition is exacerbated under limited Mek conditions where the maximum/ 
minimum values of (p are larger and smaller, respectively than the baseline conditions. Inspection 
of the dephosphorylation rates shows that an increase in dephosphorylation rate of Erkl results 
in a larger positive (p, or an Erk2p steady-state bias, and vice versa for Erk2 dephosphorylation. 
This is also reflected in the heat map representation which reveals an asymmetry, with a bias 
towards Erk2p, which is more easily activated than Erklp. This is also apparent in the time course 
(Fig. 4b (i)) where equality between Erklp and Erk2p steady state values exists for low Erkl 
dephosphorlation rate, but not for low Erk2 dephosphorylation in Fig. 4b (i). Under limited Mek 
conditions, both Erklp and Erk2p steady-state values are much smaller and activation time increases 
(Fig. 4a, b (iv)). Moreover, for limited Mek conditions, the phosphorylation/dephosphorylation 
rate parameter space has a larger region for possible competition scenarios (see largest (white) and 
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Figure 4. Competition in minimal model at baseline and limited Mek (Mek=0.2) conditions. 
Competition scenarios are studied using heat maps (middle columns) and a time course (outside 
columns) of steady state. Heat maps (ii-iii) of parameter /initial condition values are given along 
the horizontal axis (j) and vertical axis (k) are varied and the associated color at each {j, k)- 
parameter/initial condition correspond to ip = log(Erk2p/Erklp) where a positive value denotes 
more Erk2p than Erkip and a negative value indicates more Erkip than Erk2p. Solid circle and 
open circle in heat map are shown in time course of Erkip (blue) and Erk2p (red) in (i, iv) where 
a solid circle corresponds to a solid line and an open circle corresponds to a dotted line, (a) Effects 
of phosphorylation rates on Erkl/2p. Values of {k2,k&) are varied, solid circle (solid line) is high 
k2 — 0.6, low fee = 0.1 and open circle (dotted line) is low ^2 = 0.1, high fcg = 0.6. (b) Effects of 
dephospho rylation rates on Erkl/2p. Values of (feg,feio) are varied, solid circle (solid line) is high 
fcg = 0.6, low feio = 0.1 and open circle (dotted line) is low fcg = 0.1, high fcio = 0.6. (c) Effects 
of Erkl and Erk2 initial conditions. Values of total Erk concentrations (Erklr, Erk2T) are varied, 
solid circle (solid line) is high Erkly = 2, low Erk2T = 0.5, and open circle (dashed line) is low 
Erklr = 0.5, high Erk2T = 2. 



200 




200 



smallest (black) values of 99, corresponding to high Erk2p and Erkip dominance in Fig. 4a, b (iii)). 
This suggests that for certain phosphorylation/dephosphorylation rates, a limited stimulus would 
more strongly favor an Erkip response than non-limited stimuli. We investigate how the total 
amount of Erk may vary across cells affects activation states (Fig. 4c) and observe that for small 
total Mek, as well as baseline conditions, the initial amount greatly affects the steady-state value 
and it gives the expected result that as total Erkl in the system increases, the total Erkip also 
increases (Fig. 4c). Unlike the phosphorylation/dephosphorylation cases, however, the (p indicator 
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value at a given total (Erkl, Erk2) point does not change as Mek becomes limited (heat map 
indicator colors are nearly identical). 

The minimal model provides a number of insights, specifically that phosphorylation/dephosphorylation 
rates play an important role in the steady-state behavior of Erklp and Erk2p, whereas under lim- 
ited Mek conditions, the parameter rate space suggests there is a stronger effect (larger \ip\) on the 
response. The value of the initial conditions can induce a Erk2p or Erklp dominated response and 
limited Mek alters the steady-state value of this response. More generally, this model provides a 
simple framework for gaining insight into the components which control the competition between 
Erkl and Erk2 for its kinase Mek, and we provide a indicator ip for giving Erklp or Erk2p cell 
response. 

Complete system. The minimal model provides preliminary insights into the behavior of Mek 
and Erk isoforms. Even though the analysis of the minimal model is useful, further analysis of the 
complete system must be performed to acquire information about the role of translocation between 
the two cellular compartments (cytoplasm and nucleus). 

As shown in Fig. 5a, the equilibration dynamics and the steady state of Erklpn and Erk2pn in 
HeLa and NIH 3T3 cells are most sensitive to the import shuttling rates of Erklp, Erk2p, and 
M-Erk2 (54,57,53); and the only sensitive biochemical reaction rate is the dephosphorylation rate 
of Erk2pn {ki2)- The primary difference between the cell types is the sensitivity of Erk isoform 
shuttling speeds. Specifically, HeLa cells are more sensitive to Erk2p import shuttling whereas NIH 
3T3 cells are more responsive to changes in Erklp. Overall, both cell types are quantitatively most 
sensitive to the shuttling of Erk2. This analysis also highlights the similarities and differences in 
sensitivity between initial conditions of the minimal model and the complete system (Fig. 5b). The 
intermediary complex of Mek-Erk2 are sensitive in both models and the nuclear complexes are also 
sensitive. The drastic difference in steady state concentrations between the two cell conditions (lim- 
ited Mek and baseline) that was observed in the minimal model does not transpire in the complete 
model. As shown in Fig. 5a, b, the sensitivity index of the most important parameters are nearly 
100 times those of the initial conditions. Investigating the stimulus-response curve reveals that 
unlike the minimal model Mek(O) = 1 does not saturate Erk2pn and requires a stronger stimulus 
(compare Fig. 3c and Fig. 5c). The sensitivity analysis was repeated for limited Mek conditions, 
which did not present any difference in the importance of parameters or initial conditions. 

Given the importance of the shuttling rates on Erklpn and Erk2pn concentrations, we investigate 
the response of HeLa and NIH 3T3 cells to different conditions affecting shuttling. For both cell 
types, we find that increasing the Erklp shuttling rates (Fig. 6a), Erk2p shuttling rates (Fig. 6b), 
and Mek shuttling (Fig. 6c) by an order of magnitude increases the steady-state and the activation 
time. While the sensitivity analysis suggests that NIH 3T3 cells are more sensitive to Erklp 
shuttling, the differences between Fig. 6a (i-ii) are slight. In NIH 3T3 cells, the activation time of 
Erklpn is shorter and the steady state is slightly larger. As shown in Fig. 6a, b, Erklp shuttling rates 
only affect Erklpn shuttling and Erk2p shuttling only affect Erk2pn steady state. By increasing 
and decreasing the Erk2p shuttling rates, we find a quantitatively larger difference in steady state 
value of Erk2j»n. This difference confirms why Erk2p shuttling is most important in the complete 
system. By examining the effects of Mek shuttling on Erkl/2pn steady states, it is clear that Erk2pn 
responds more than Erklpn both in activation time and steady state value. While it is well known 
that Mek has a nuclear exclusion sequence that exports it from the nucleus, studies have suggested 
that Mek may play an active role in transport of Erk |40| HI] . Our findings that the shuttling of 
Mek is less sensitive than Erklp/Erk2p, yet result in the same effect on steady state and activation 
time (compare Fig. 6b and Fig. 6c) is somewhat counterintuitive. Evidence suggests that the 
Mek mediated contribution to Erk export is at most marginal under physiological conditions, only 
contributing to the steady-state in unphosphorylated conditions [8J. This may be due to the very 
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low rate supported by the CRMl mediated export. Despite the numerous studies to understand 
the role Mek plays in regulating Erk, our analysis gives a juxtaposition: Mek shuttling sensitivity 
is relatively small compared to Erklp/Erk2p shuttling sensitivity while simulations predict that 
Mek and Erk shuttling have the same effect on activated nuclear Erk. 

The effects of dephosphorylation on phospho-nuclear Erk steady-states is more pronounced than 
the effect of phosphorylation rates. By varying the dephosphorylation rates of Erk2 (fci2 = 
1.4 X 10~^ baseline) a decrease by an order of magnitude remarkably affects the steady state 
value of Erk2pn. Increasing the dephosphorylation rate more than four-fold {k^ = 6 X 10~^ ) 
decreases the steady state approximately half of the baseline; however the potential for increasing 
the activation time and steady state value by changing the dephosphorylation confirm the findings 
from the minimal model. By studying the effects of phosphorylation rates on the steady-state val- 
ues, the phosphorylation rates do not exhibit a notable change in either activation or steady-state, 
which is a significant difference between the minimal and complete models. 

The complete system is consistent with the main finding that the shuttling rates of phosphory- 
latcd Erk affect the steady-state values of Erklp„ and Erk2p„ (Fig. 5a). The sensitivity analysis 
for HeLa cells suggests it is more sensitive to Erk2p shuttling whereas NIH 3T3 cells are more 
sensitive to Erklp shuttling; however our simulations do not show remarkable differences between 
the two cell types. We also observe that both cell types are responsive to Mek shuttling speeds 
although Mek shuttling parameters are not overly sensitive in analysis, which may explain the 
ongoing investigations of the role of Mek shuttling on Erk. 

Discussion. Our analysis suggests that transport of Erk is most important for determining ac- 
tivity levels. Competition between Erkl and Erk2 for Mek is induced as biochemical rates and 
total concentrations of Erk and Mek change; our minimal model and analysis may also offer an 
explanation of why the Mek-Erk stochiometry is close to one. Finally, analysis of the complete 
model highlights the unexpected result that as Mek shuttling rates are less sensitive in analysis, 
the simulation of the solution curves provides the same effects as varying the Erk2p shuttling rate. 

The differences between cell types are also illuminating: since the shuttling rates are the most 
sensitive parameter one may be led to think that differences in trafficking are at the basis for 
heterogeneity. There are two additional points that we feel are pertinent in this context: (1) the 
composition and amount of nucleoporins changes during the cell cycle, which bring about changes 
in trafficking rates and in levels of Erk activation; (2) in neuronal cells, where trafficking rates 
are maximally different between Erkl and Erk2 (unpublished data), we may expect enhanced 
difference between Erkl and Erk2 activation levels. This is notable (it is certainly enthralling 
from a neurobiological point of view) since the only clear phenotype of the Erkl knockout is 
neurophysiological. Although beyond the scope of this study, considering changes in transient 
changes in Erk isoform concentration as well as exploring the effects of downstream substrates may 
pin down the physiological implication of the shuttling mechanisms discussed here. 

In response to a given stimulus the different activation times are mostly due, not to differences 
in the phospho-dephospho rates, but to differences in shuttling speed. The FRAP experiments 
show that there is a large variability in trafficking speed, likely due to the variable density of 
nucleoporins, and this could go a long way towards justifying the different kinetics. In view of 
the large individual differences between cells, this seems to be interesting and hitherto neglected 
aspects of studies of cell-to-cell variability. Through the development of these abstractions and 
sensitivity analysis as performed might help to understand what are the most likely parameters 
responsible for the heterogeneous responses within a cell population. 

More generally, our results highlight the important aspect of spatial localization of specific iso- 
forms that has hitherto often been neglected in models and experiments. Similarly to MAPK stud- 
ies, NEAT has been found to have different behavior based on stimuli duration (sustained, pulsed 
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Figure 5. Analysis of complete system. Sensitivity coefficients for parameters (a) and initial 
conditions (b). For both HeLa and NIH cells only a limited number of parameters have a significant 
influence on the behavior of Erklpn and Erk2pn: S7, S4, se, S7t-, fci2, S2, S3- Sensitivity to initial 
conditions reveals that for both cell types all initial conditions have a comparable impact on the 
observed dynamics. The most sensitive initial conditions, however, are much less sensitive than 
most sensitive parameters, (c) Effects of phosphorylated Mek on output concentrations Erklpn and 
Erk2pn. 



etc); notably, NFAT isoforms in the Calcium/NFAT pathway that regulate T cells have different 
timed responses which may be explained by tissue specific subcellular localization (Nir Friedman, 
personal communications). These are not isolated cases: even Fox-1 homologs are tissue-dependent 
and a theoretical study of their nuclear localization, such as this one, may determine how these 
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Figure 6. EfTects of Erk and Mek shuttling, and dephosphorylation rates for complete model at 
baseline parameter values and initial conditions. Solid line denotes HeLa cell baseline and dotted 
line denotes NIH 3T3 cell baseline. a,b) Effect of Erklp and Erk2p shuttling speeds. The shuttling 
parameters, S4 and S7, are varied by an order of magnitude from their baseline values. The solution 
curves for varying Erklp shuttling rates is given by slow (dots) and fast (cross) whereas Erk2p 
shuttling rates is given by slow (squares) and fast (diamonds) can be compared to their respective 
baseline values, c) Effects of Mek shuttling rates. The Mek import shuttling parameters, S2, S3 and 
Se, are varied from their baseline values. The time course for slow Mek shuttling is denoted by stars, 
fast Mek shuttling is represented by upside-down diamonds lines, d) Effects of dephosphorylation 
rate of Erk2p. The phosphorylation parameter, fci2 is varied from its baseline values. The solution 
curves for slow rates (fci2 = 0.014) are denoted by circles, fast rate (fci2 = 0.6) are represented by 
triangles. 



isoforms regulate activities of neural tissue-specific splicing. Thus we feel that differential shuttling 
of signaling molecule isoforms coupled to competition between them may have a general role on 
cellular decision making processes. 
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Supplementary Material 



There are two sections which include all model equations, minimal model parameters and initial 
conditions and Fig. SI. 



Models 

Equations. The biochemical reactions of the full system given in Table 1 are translated into a 
system of ODEs following mass action kinetics: 



d [Erkl 
d [Mek 

d [M-Erkl 
d [Erklp; 
d [Erk2 

d [M-Erk2 
d [Erk2p; 
d [Erkl„ 
d [Mek„ 

d [M-Erkl„ 
d [Erklpn 
d [Erk2„ 

d [M-Erk2„ 
d [Erk2p„ 



/dt = -fci [Erkl] [Mek] + /c_i [M-Erkl] + fcg [Erklp] + s_i [Erkl„] - si [Erkl] , 
/dt = -ki [Erkl] [Mek] - fcg [Erk2] [Mek] + (fc_i + fca) [M-Erkl] 

+ (fc-5 + ke) [M-Erk2] + s_2 [Mek„] - sa [Mekl] , 
/dt = ki [Erkl] [Mek] - (A:_i + ita) [M-Erkl] + s^g [M-Erkl„] - S3 [M-Erkl] , 
/dt = k2 [M-Erkl] - kg [Erklp] + s_4 [Erklp„] - S4 [Erklp] , 
/dt = -fcg [Erk2] [Mek] + k-5 [M-Erk2] + fcn [Erk2p] + S-5 [Erk2„] - sg [Erk2] , 
/dt = its [Erk2] [Mek] - {k-5 + ke) [M-Erk2] + s_6 [M-Erk2„] - se [M-Erk2] , 
/dt = ka [M-Erk2] - fen [Erk2p] -|- S-7 [Erk2p„] - st [Erk2p] , 

/dt = -fes [Erkl„] [Mek„] + k-3 [M-Erkl„] + feio [Erklp„] + si [Erkl] - s_i [Erkl„] 
/dt = -ks [Erkl„] [Mek„] - fcg [Erk2„] [Mek„] + (fc_3 + ki) [M-Erkl„] 

+ (fc_7 + kg) [M-Erk2„] + S2 [Mek] - s_2 [Mek„] , 
/dt = fes [Erkl„] [Mek„] - (fc_3 + kt) [M-Erkl„] + S3 [M-Erkl] - s_3 [M-Erkl„] , 
/dt = fe4 [M-Erkl„] - feio [Erklp„] + S4 [Erklp] - s_4 [Erklp„] , 
/dt = -fer [Erk2„] [Mek„] -I- k-7 [M-Erk2„] -f- fei2 [Erk2p„] + ss [Erk2] - s_5 [Erk2„] 
/dt = fer [Erk2„] [Mck„] - (k-7 + fes) [M-Erk2„] -I- ss [M-Erk2] - s_6 [M-Erk2„] , 
/dt = ks [M-Erk2„] - fei2 [Erk2p„] + S7 [Erk2p] - s_7 [Erk2p„] . 



Note, continuous shuttling is dependent on Mek activation [8]; therefore, we assume that Mek 
activation is sustained in the system such that total active Mek is constant, 

Mekr = Mek + Mek„ + M-Erkl + M-Erkl„ + M-Erk2 -|- M-Erk2„. 

Since this study only focusses on the short timescales of Erkl and Erk2 activation and not ensuing 
gene expression processes, we assume that total Erkl and Erk2 are also constant and satisfy the 
conservation relations, 

ErklT = [Erkl] + [Erkl„] + [Erklp] + [Erklp„] + [M-Erkl] + [M-Erkl„] , 
Erk2T = [Erk2] + [Erk2„] + [Erk2p] + [Eiklpn] + [M-Erk2] + [M-Erk2„] , 

where Erklj- and Erk2'p are total Erkl and Erk2, respectively. 

The reactions of the minimal model are given in Table 1 and the relevant reactions are identified 
by only cytoplasmic reactions, corresponding to the first seven differential equations in the full 
model system, except the terms with a coefficient of Sj are disregarded. Thus the model in Eqs. 1 — 7 
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become: 

(2) d [Erkl] /dt = -fci [Erkl] [Mek] + fc_i [M-Erkl] + fcg [Erklp] , 

(3) d [Mek] /dt = -fci [Erkl] [Mck] - fcs [Erk2] [Mek] 

(4) + (fc-i + ki) [M-Erkl] + + fee) [M-Erk2] , 

(5) d [M-Erkl] /dt = fei [Erkl] [Mek] - (fe_i + fea) [M-Erkl] , 

(6) d [Erklp] /dt = fea [M-Erkl] - feg [Erklp] , 

(7) d [Erk2] /dt = -fcs [Erk2] [Mek] + fes [M-Erk2] + fen [Erk2p] , 

(8) d [M-Erk2] /dt = fes [Erk2] [Mek] - (fcs + fee) [M-Erk2] , 

(9) d [Erk2j3] /dt = fee [M-Erk2] - fen [Erk2p] , 

and conservation relations: 

Mekr = [M-Erkl] + [M-Erk2] + [Mek] , 

Erklr = [Erkl] + [Erklp] + [M-Erkl] , 

Erk2r = [Erk2] + [Erk2p] + [M-Erk2] , 

where Mek^, Erkly and Erk2T G M are total Mek, Erkl and Erk2, respectively. We introduce the 

normalised variables [Ei^l] = [Erkl]/ErklT, [Erklp] = [Erklp] /Erklr, [^rk2] = [Erk2]/Erk2T, and 

[Erk2p] = [Erk2p]/Erk2T. Using the conservation relations above and the new normalised variables 
we obtain the following equations: 



d 



d 



d 



Erkl 

Erklp 
Erk2 
Erk2p 



where 

Mi = l- 



/dt = 
/dt = 
/dt = 
/dt = 



-fciErklT 

k2Mi - kg 
-A;5Erk2T 
k6M2 - kii 

Erklp 



Erkl 

Erklp 
Erk2 
Erk2p 



Erkli 



/ Mckr 
V Erk2r 



Ml - nM2 + k-iMi + kg 



Erklp 



n~^Mi - X2 ) + k5M2 + fell 



Erk2p 



M2 = l- 



Erk2 



Erk2p 



and 



7^ = 



Erk27 



ErklT 

and the system has been reduced by three equations. For simplicity, in the main section of the 
article, we drop the tilde notation from the variables. 

Estimation of transport parameters and initial conditions. Before photobleaching, the cells 
are assumed to be in equilibrium and that influx and efflux of the Erk species are equal. After 
bleaching, the flux continues in and out of the cell, and fluorescence is measured. We introduce two 
equations to describe the bleached nuclear Erk species (F^) and recovered fluorescence of nuclear 
Erk species (-Fn): 



dt 



and 



dt 



where s-j is the rate at which bleached species leaves the nucleus and a is the rate at which 



unbleached, fluorescent species enters the nucleus. At equilibrium 



dt 







a 



[Fn] 

a/s-j ' 

Then the solution is 

Fn = 1- e~*'*-J'. If at t = T, and F* is the measured fluorescence (F„ = F*), then F* = I — e""^'*-^ 
and the nuclear export rate is given by s-j = — — -. Here, the F* is the recovered fluorescence 
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we have [F„] = a/s-j. As before, we normalize by the steady-state of Fn (a/s-j): F„ 
giving [Fji] = -^Fn, so at equilibrium F„ = 1, 



and ^ 



S-i 



S-jFn 



at the time constant r. In NIH 3T3 cells we have for starved Erkl r = 653s, stimulated Erklp 
r = 266s, starved Erk2 r = 178s, and stimulated Erk2p r = 84s when F* « 0.75 normalized 
fluorescence [5] . Given that the sizes of HeLa and NIH 3T3 cells are similar (the ratio of molecules 
in cytoplasm/nucleus is 1.5 it 0.2 for Erk and 1.1 it 0.72 for Erkp) the import rates (sj) can be 
estimated as export rate divided by the molecular number ratio, as calculated by [29]. Without the 
shuttling rates of Mek for NIH 3T3 cells, we assume the values as calculated by [2^ for HeLa cells. 

Initial conditions of the transport model are the same as the full model. The initial conditions of 
the biochemical model follow the normalization such that Erkl(0)=l, Erk2(0)=l and Mek(0)=l, 
and the total concentrations of Erkl (Erkly) and Erk2 (Erk2r) are set to maintain that Erk2 is 
four times more abundant than Erkl. 



Total Nuclear Erk 

Supplementary Fig. 1. Time course of total nuclear Erk for HeLa and NIH 3T3 cells A. Full model, 
(i) Time course of total nuclear Erkl (Erklp„ + Erkl„) and Erk2 (Erk2p„ + Erk2„) at baseline 
conditions, (ii) Time course of total nuclear Erkl (Erklp„ + Erkl„) and Erk2p (Erk2p„ + Erk2„) 
at equal Erkl/2 (Erkl = Erk2 = 0.5). (in) Time course of total nuclear Erkl (Erklp^ -|- Erkl^i) 
and Erk2 (Erk2p„ + Erk2„) at Limited Mek (Mek = 0.2). B) Transport model. Transport model 
was simulated at baseline parameter values for HeLa cells (solid) and NIH 3T3 cells (dashed). 
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